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LOCAL ERROR CONTROL AND ITS EFFECTS 


ON THE OPTIMIZATION OF ORBITAL INTEGRATION 


INTRODUCTION 

Many computing systems which are used to determine the orbits of artificial 
earth satellites require numerical integration. Due to advances in the theory of 
perturbations and the concomitant increase in the complexity of the mathemati- 
cal model, highly efficient integration techniques are desirable and, in certain 
cases, necessary. 

The orbital equations under consideration are of the form 


x 


-fix 

IIXII 3 


+ p 



( 1 ) 


in three space variables, where nxn = (x 2 + y 2 + z 2 )^, P is the perturbation 
function and ^ is a constant. We will assume throughout that p is fairly com- 
plicated so that the efficiency of the integration is proportional to the total num- 
ber of derivative evaluations. (Also, for simplicity we assume that P is inde- 
pendent of the first derivative x, i.e., P = P(t, x) ; see Appendix B.) 

The numerical method under consideration is a multistep process, i.e., a 
method which approximates the solution of (1) by using formulas of the form 

k k 

/ a.x = h 2 

/ . 1 n 1 

1=0 i =0 


where h is the stepsize, x. = x (t Q + ih), and a. , /3. are constants. Formulas 
of this type define a linear k-step method (see Henrici 1962, pp. 295 ff.) , and 
are usually derived by the integration of a polynomial approximation of the 
derivative x. 
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/3.x 

' l n - i 


n - k, k + 1 . 


( 2 ) 
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Associated with (2) is a local truncation error of the form 


R n = R(t n ) = C h p + 2 x<p + 2 >(^) , 


( 3 ) 


where ^e[t n _ k , t n ], C is a constant, and p is an integer called the order of the 
method, all depending on k , the number of "backpoints" used. 

It has been suggested (Soar, 1964) that controlling the local error by varying 
the order and stepsize during the integration of the equations of motion (1) may 
yield gains in efficiency without sacrificing accuracy. For example, it is clear 
that the magnitude of (3) is sensitive to variations in p or h . Suppose that 
during the integration, the magnitude of R n becomes insignificant relative to 
the calculations being performed. Then an increase in the magnitude of h or a 
decrease in the magnitude of p, either separately or in combination, may in- 
crease the efficiency of the integration without sacrificing accuracy. 

The purpose of this study then, is to develop techniques to automatically 
control the magnitude of R n during the integration by varying the parameters 
p and h , and to examine the effects* of these controls on the efficiency and 
accuracy of the process. 

Before discussing methods of estimating and controlling the local error, we 
formulate a commonly used integration model which was used to obtain the 
results. 


THE INTEGRATION FORMULAE 

We begin by remarking that in equation (2) , if /3 Q / 0 , then knowledge of the 
solution x n is required on both sides of the equation and is, in general, not 
explicitly solvable. Equations of this type (closed form) however, have smaller 
associated truncation errors as well as desirable stabilizing characteristics; 
(Hull and Creemer, 1963). The well-known predictor-corrector algorithm 
utilizes formulas of this type by first computing an initial (predicted) approxi- 
mation of the solution using a similar formula with /3 Q = 0, then using the closed 
form iteratively until convergence is achieved. It can be shown that for suffi- 
ciently small h , the successive "corrected" values obtained by this process 


*It should be noted that local error control does not in general yield a quantitative appraisal of the 
accumulated error. A qualitative control of the accumulated error is possible though, and in most 
cases is sufficient. For a discussion on accumulated error estimations, see Henrici, 1962. 
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converge to the unique solution of the closed form equation provided the function 
being integrated is sufficiently smooth. 

Consider now the predictor-corrector formulas which are derivable from 
Newton's backward difference interpolation polynomial (Henrici 1962, pp. 
291-293): 


(Pred) x + 2x , + x 

v ' n n ~ 1 n - . 


V 2 x = h 2 

n 


! + 12 y2 + 12 y3 + 240 V4 


+ 40 VS + 


+ cr Vi 

q 


n - 1 


(4) 


(Corr) x n - 2x n _j + x n _ 2 - V 2 x n = h 2 

1 


1 - V + T2 V2 - 2*5 V4 


2 4 o V s + • • • + a* V«» 


(5) 


These are the Stormer- Cowell formulae expressed in terms of backward 
differences. The local truncation errors associated with these formulae are 
given by (3) where p = q + 1 and C : ^ qtl orcr* tl , We note that if the param- 
eter q is fixed, (and hence the order p) , then (4) and (5) can be written in the 
form of (2) by expressing the backward differences in terms of the ordinates x.. 

In particular, for q = 5 we have the Stormer-Cowell formulas of order 6: 


V 2 , 


,2 


240 


317 x - 266 x 0 + 374 x - 276 x 

n-1 n~2 n — 3 n - 4 


+ 109 X - “ 3 X , 

n "* 5 n — 6 


(4)’ 


V 2 , 


,2 ^ 


240 


18 x + 209 x , + 4 x + 14 x ,-6x . + x . 

n n-1 n ~ 2 n - 3 n~4 n~5 


(5)' 


It is clear that the coefficients in (4)' and (5)' depend on the choice of order, 
so that varying the order during the integration would mean producing a new set 
of coefficients for each p. It is for this reason that formulas of the type (4) and 
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(5) , where variations of the order can be made simply by varying the number 
of terms retained, were used in this study. (*) 

We note that before any of the above formulas could be used, a set of 
"starting" values of the solution must be computed. (**) For example, for 
equations (4)' and (5)', if the values x. (and hence x.), i = 0, 1, •••5 are 
known, then (4)' could be used to obtain a "predicted" value x£ (n =6) and (5)’ 
to obtain the successive corrections x^ j , j = 1, 2,3, • • • until convergence to 
same criterion 8 is achieved, i.e., until 


c. . , c. . 

V ■ V I - 8 • 


The process could then be repeated with n = 7, 8, 


LOCAL ERROR ESTIMATION 

Any considerations concerning the control of the local error depend on the 
capability of obtaining a reasonable estimate of (3) during the integration. One 
of the most widely used approximations is based on the fact that if f(x) is an 
n-times continuously differentiable function, then there exists a 8 , 0 < 6 < 1, 
such that 


A n f ( x ) _ d n f 

(Ax/ 1 dx n 


(x + ni 9 (Ax)), i.e., 


if Ax is sufficiently small, we can approximate a high order derivative by a high 
order difference. In particular, we can write 


R = R 


(t n ) = ChP + 2 x<p + 2 >(£) = Ch p + 2 i< p >(£) = Ch 2 V p x n (6) 


(*) In actual practice, a modification of equations (4) and (5) known as the "summed” form of the 
integration formulae was used. This modification is formulated in Appendix A. 

(**) Experimentation has indicated that recently developed high order Runge-Kutta type formulas 
arc particularly suited for such purposes; see Shanks, 1966. 
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If suitable predictor-corrector formulas are used, another technique 
available is based on comparing the predicted value with the finally accepted 
"corrected" value x^; in particular, it can be proven (Henrici, 1962) that for 
equations (4) and (5) , 


R 

n 



(7) 


where K is a constant. This technique is known as Milne's method of estimating 
the local error. 


LOCAL ERROR CONTROL 

We consider now a specific technique which would enable one to control the 
magnitude of the local error during the integration by varying the parameters 
p and h. Let Tj and T 2 be tolerances specifying the desired upper and lower 
bounds on the local error, so that for any n, the local error must satisfy 

T 2 < | R n 1 £ Tj, where T 2 < Tj . (8) 


Controlling the local error by varying the order can then be accomplished 
by determining whether condition (8) is satisfied for each n, i.e., if for some n, 
R n < T 2 decrease the order, or if R n > Tj, increase the order. 

We note that varying the order alone is generally not sufficient for an effec- 
tive control. For example, if for some n, R n > Tj and h is large, it may not 
be possible to satisfy (8) with any p, in which case the stepsize must be de- 
creased. (Also, the danger of numerical instability increases considerably with 
large p; see references 5 and 6.) On the other hand, R n < T 2 and p small would 
indicate that a larger stepsize could be used. 

Controlling the local error by varying the stepsize, in addition to varying 
the order, can be accomplished as follows: Let Lj and L 2 be limits specifying 
the desired upper and lower bounds on the order, i.e., 

L 2 < p < Lj where L 2 < Lj ; (9) 

then, if at some point during the integration, p > Lj , decrease the stepsize, or 
if P < L 2 , increase the stepsize. 
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A local error control designed in this way has the property that the param- 
eters Tj, T 2 , Lj, L 2 completely govern the degree and type of control; i.e., 
the parameters Tj and T 2 can be selected so as to effect any degree of control 
from none, to a continuous step by step control; likewise, the parameters L x and 
L 2 can be selected so that the control involves varying the order alone, varying 
the step alone, or varying both. For example, if Lj = L 2 , the control depends 
solely on variations in the stepsize. 

We remark that unlike varying the order, changes in the magnitude of the 
stepsize during the integration are usually not easily accomplished since a 
"memory" of equally spaced points is required at every step during the inte- 
gration. For this reason a common technique is the "halving-doubling" method 
where an "increase" or "decrease" in the stepsize is either half, or double the 
current stepsize. Then increasing the stepsize presents no problem if a suffi- 
cient number of backpoints are retained. Decreasing is usually accomplished by 
some interpolation technique or by using a single step method. 


OPTIMIZATION OF STEPSIZE 

Restricting the variations of stepsize, h , by a constant factor does not, in 
general, yield optimal stepsizes*, and hence the initial choice of interval could 
have a substantial effect on the total number of integration steps. For example, 
suppose the halving-doubling method is being used and the order is held fixed. 

If for some n , R n > T x , the stepsize must be decreased. Let h op( .be the opti- 
mum (largest) value of h for which (8) is satisfied. We could then have the 
situation h/2 < h < h. Because of our restriction h/2 would be used, although 

opt 

this would result in more integration steps than would be required if the optimum 
stepsize were used. Hence a variation of h which would better approximate its 
optimum value is desirable. One technique available (see Reference 4 for details 
and applications in the case of single-step methods) is the computation of the 
stepsize using the local error estimate, i.e., let a be the "allowable" local 
error for each step , where T 2 < a < T r Variations in h can then be computed 
by using the relationship between a and R n = Ch p + 2 x (p + 2 >(£) , viz., 



a 

1/p+2 ^ 

'crh p + 2 ~ 


_Cx( p + 2 >(£)_ 


R 

L n J 


(10) 


*The largest stepsize which allows a prescribed local error at a given point. 
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so that if cr < R n , the stepsize is decreased, or if R n < a the stepsize is in- 
creased, where the variation is approximately optimal with respect to the choice 
of cr . 


OPTIMIZATION OF ORDER 

As in the case of varying the stepsize, we see that varying the order by 
constant factor need not yield the optimal* order. For example, suppose a 
variable order, constant step control is being used and the order is varied by 
±1. If for some n, |R I < T, , the order must be decreased. Let p be the 

1 n 1 2 K op t 

optimum (smallest) value of p for which (8) is satisfied. We could then have the 
situation 


P op t < P-1 < P- 

Because of our restriction p - 1 would be used, although this would result 
in more calculations than would be required at the optimum order. Also, if we 
were using a variable step control in addition to varying the order , it may occur 
that 


Popt < L 2 < P - 1 , 

in which case the stepsize would be increased if p opt were used. One method 
which could be used to obtain the optimum order is to test the local error for 
various .orders until the optimum is established. For example, if the local error 
estimate used is given by (6) , then one could test Ch 2 VPx n for various p until 
the smallest value satisfying (8) is found. 


GENERAL EFFICIENCY CONSIDERATIONS 

Let us examine some possible effects, advantages and disadvantages of the 
above mentioned controls on the efficiency of any integration, which will 
generally depend on the following: 

(a) Minimizing the number of integration steps, 

(b) Minimizing the number of corrector iterations required at each step 


* 


The smallest order which allows a prescribed local error. 
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(c) Minimizing other computational efforts required by the integration 

formulae being used, such as retaining only the significant terms in (4) 
or (5) during the integration. 


Further, we must consider 

(d) The computational effort and time involved in controlling the local error, 
such as producing the required "memory" when changing the stepsize. 


We now examine the following controls: 


(I) A variable order, fixed step control 

Since varying the order involves only varying the number of terms retained, 
this control has little effect on (d). Equation (7) however, expressing the re- 
lationship between the local error and the "predicted-corrected” difference, 
indicates that controlling the local error may substantially effect the total num- 
ber of corrector iterations. In particular, an increase in order at some point 
during the integration where the local error is increasing could minimize any 
increase in the number of required iterations. 


Another possible advantage of this control is its obvious effect on (c), and 
that given a stepsize, the order required to satisfy a given criterion on the local 
error is automatically selected. 

A clear disadvantage of this control is that it has no effect on minimizing 
the number of integration steps. 


(II) A fixed order, variable step control 


This control affects the number of corrector iterations for the same reason 
as (I). It is clear that (a) and (d) are affected, and we must consider the "trade- 
off" between these two, i.e., we must consider the gain in efficiency due to the 
larger stepsizes against the cost in changing the stepsize. Furthermore, we 
have discussed two possible methods for varying the stepsize, viz., (i) halving- 
doubling, and (ii) step computation using the local error estimate. In both cases 
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the necessary "memory" could be obtained by interpolation, but in method (ii) 
all the backpoints must be computed both when decreasing and increasing, as 
opposed half, or none, for method (i). On the other hand, we have seen that 
method (ii) could be more effective in minimizing the integration step than 
method (i). 


Besides its possible effects on (a), another advantage of this control is 
that given an order, the stepsize (possibly optimum) required to satisfy a given 
criterion on the local error is automatically selected. 


(Ill) A variable order, variable step control 


The control combines the effects on (I) and (II). A possible disadvantage 
of this control is that in allowing the order to vary before changing the step, 
smaller stepsizes may result, thereby adversely affecting (a), (as opposed to 
control (II)). 

A clear advantage of this control is that both the stepsize and order required 
to satisfy a given local error criterion are automatically selected. 


EFFECTIVE ERROR CONTROL IN ORBITAL INTEGRATION 


In general, the effectiveness of a local error control during the integration 
of a particular orbit depends on the degree and rate of change of the derivative 
x*(p + 2 )(£) t (and hence the local error), during a revolution for a fixed p and h . 
This change is governed by the orbital parameters a, the semi-major axis, 
and e, the eccentricity (Soar, 1964). In particular, for orbits with eccentri- 
city near zero, it can be expected that the local error variation will be 
small over a revolution, and that as the eccentricity becomes bounded away 
from zero, this variation becomes more pronounced. Moreover, the larger the 
semi-major axis is for a particular satellite, the slower the rate of change 
of the derivative. 


Assume now that during the computation of a particular orbit, some local 
error criterion is to be satisfied over the entire range of integration, (which may 
involve many revolutions), i.e., we wish to bound the local error from above, so 
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as to restrict, at least qualitatively, the propagation of truncation error. As- 
sume also that the integration method used is of some fixed order, (as is 
generally the case when formulas of the type (4)' or (5)' are used) , and that the 
stepsize is to be specified. In such a situation, one could insure the error cri- 
terion will be satisfied and yet obtain some degree of efficiency, by integrating 
over a single revolution with various stepsizes, and selecting the largest step- 
size which satisfies the criterion over the entire revolution. For orbits with 
e = 0, this process could yield the optimal stepsize for the given order. On the 
other hand, for orbits with e > e > 0, the stepsize selected in this manner would 
correspond to bounding the local error at its maximum value over the revolution 
(e.g., at perigee) , and could result in needless computation for those portions 
of the revolution where the local error is smaller, or at its minimum, (e.g., 
near apogee) . 

Consider now the effects of a continuous variable step control during the 
integration. In both cases cited above, i.e., e = 0, ore > e > 0, an initial 
stepsize (possible optimum) satisfying the local error criterion would be com- 
puted automatically, and furthermore, the stepsize would vary according to 
variations in the derivative. It will be seen in the numerical results that in both 
cases , under suitable conditions , such a control has a considerable effect on the 
efficiency. In particular, it will be seen that if a sufficiently large range 
(Tj , T 2 ) in (8) is used, so that the local error is allowed to range between these 
two limits before any interval modification is performed, the gains in efficiency 
due to the larger stepsizes will overshadow any loss due to changing the stepsize. 

If, in addition to the above, the order to be used could be specified, (as is 
generally the case when formulas of the type (4) or (5) are used) , then it may 
be possible to obtain an optimal stepsize (or stepsizes) - order combination by 
integrating over a single revolution with various orders and letting the stepsize 
vary during the revolution. In particular, one could find that order which, 
together with a continuous step modification, effects the most efficient inte- 
gration. Examples of such a procedure will be given in the numerical results. 

Finally, let us examine the effects of a continuous variable order control. 
We note again that since we are assuming that the number of evaluations of the 
derivative governs the overall efficiency , the effects on this control or (c) (in 
the previous section) , and its possible adverse effects on (a) when used in con- 
junction with a variable step control, make this control ineffective from this 
point of view. If, however, a fixed step method is used, it will be seen that the 
effect of this control on the number of corrector iterations can result in con- 
siderable gains in efficiency. 
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NUMERICAL RESULTS 

In this section the results obtained by applying the various local error con- 
trols during the integration of three selected orbit types will be examined. The 
formulation of the actual equations of motion (1) used is presented in Appendix B, 
and a description of the computer program can be found in Reference 9. 

Before proceeding, we make the following remarks concerning the numerical 
computations: 

(i) In actual practice, for simplicity of computation, instead of using an ap- 
proximation for Rn (given by (6) in controlling the error, a bound on the local 
error given by 


Un = Ch 2 V k " 2 x 

n 


was used, where V k ' 2 is the last backward difference retained in the computa- 
tions; (see Appendix A). Since p = k + 1 , we have 


Rn < Un 


and thus the qualitative results obtained by using Un are the same as if Rn were 
used, and the same quantitative results could be obtained by using Rn and a 
smaller T r 

(ii) The number of derivative evaluations for each integration was obtained 
by multiplying the total number of steps taken by the average number of predictor- 
corrector iterations. 

(iii) The error estimates tabulated were obtained by integrating (1) with 
P = 0, and comparing the results with the Keplerian solution. Since the effect 
of the perturbation function P on the accumulation of error is generally small, 
the error estimate obtained in this manner can be assumed to be a good estimate 
of the actual error generated. 

(iv) All computations were performed on the Uni vac 1108 computer using 
double-precision arithmetic. In the computations the following units were used: 

unit of length = 6378.388 Km. 
unit of time = 13.4472 Mins. 

We begin by tabulating (Table I) the results obtained by integrating the 
orbits with formulas of orders 7, 9, 11 and 13. These integrations were carried 
out with various stepsizes, where in each case the stepsize used was increased 
until Un (see remark (i)) failed to satisfy the inequality 

| Un | < T t (12) 


11 



over the entire range of integration, which was chosen arbitrarily to be 4000 mins. 
We have tabulated only the last 2 stepsizes tested, indicating the largest stepsize 
which passed the criterion, and the first stepsize which failed. 

Table I 

Fixed Order - Fixed Step 


Tj = .5 x 10 8 S = .1 x 10 10 (predictor-corrector tolerance) 


a 

e 

Order 

Stepsize 

(Mins.) 

No. of 
Der. Eval. 

Estimated 

Error 

6.7 

.003 

7 

5.0 

798 

.4 x 10” 9 




7.0* 

569 

.4 x 10" 8 



9 

15.0 

262 

.3 x 10" 8 




17.0* 

231 

.8 x 10" 8 



11 

24.0 

160 

.9 x 10" 9 




26.0* 

147 

.2 x 10" 8 



13 

22.0 

173 

.3 x 10" 11 




24.0* 

158 

.9 x 10" 10 

1.15 

.075 

7 

.4 

9998 

.8 x 10" 7 




.5* 

7998 

.4 x 10" 6 



9 

.9 

4440 

.2 x 10" 6 




1.0* 

3996 

.5 x 10" 6 



11 

1.2 

3327 

.5 x 10" 7 




1.3* 

3070 

.1 x if)" 6 



13 

1.5 

3081 

.1 x 10~ 8 




1.6* 

3043 

.1 x 10" 7 

8.5 

.878 

7 

.30 

13344 

.1 x 10" 6 




.35* 

11483 

.4 x 10' 6 



9 

.40 

10005 

.4 x 10" 7 




.45* 

8902 

.1 x 10' 6 



11 

.30 

13340 

.9 x 10" 10 




.35* 

11433 

.4 x 10" 9 



13 

.20 

19992 

.4 x 10" 11 




.25* 

15992 

.3 x IQ" 11 


♦Integration with this stepsize failed the criterion (12) for some n. 
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Concerning the results in Table I, we make the following observations: 

(i) For a given local error criterion, increasing the order does not gen- 
erally imply a larger "allowable” stepsize. Since the propagation of error 
(both truncation and round-off) influences the "smoothness" of the higher order 
differences (and hence our local error estimate), one would expect this behavior 
both from the inaccuracies in the differences as well as an unstable p and h com- 
bination. We also note an increase in accuracy with the higher orders. A table 
demonstrating this behavior (for the case e = 0) can be found in Reference 6. 

(ii) The stepsizes resulting from our local error criterion during the in- 
tegrations correspond to bounding the local error at its maximum (perigee), 
although the local error for most of the integration, particularly for the case 
e = .87, was considerably smaller than our upper bound. 

In Table II the results obtained by integrating our test orbits with a variable 
step control are tabulated. In particular, the stepsize was varied during each 
integration forcing Un to satisfy 

Tj < | Un | < T 2 (13) 

for all n. The stepsize modification techniques used were the halving-doubling 
and optimum step computation where the "allowable" local error used is designated 
by o- (see table below); i.e., the step computation is given by 

h o P t = [ah p+ 2 /Un] 1 / P +2 . 

Except for those cases indicated by an asterisk, all integrations were performed 
with an initial stepsize (chosen arbitrarily) of .42 mins. = 1/2 5 in internal units. 

We have indicated the stepsize or stepsize range which occurred during the 
multistep integration as a result of our local error control. 

Comparing the corresponding results in Tables I and II, we make the follow- 
ing observations: 

(i) For the "low e" cases, (e = .003, .075), the stepsizes selected as a re- 
sult of the valuable step control are comparable to the "optimum" stepsizes 
found by trial and error in Table I. 

(ii) The dependence on the initial choice of stepsize in the halving- 
doubling type of step modification is demonstrated by the cases in which the 
initial step was selected as a result of foreknowledge of the "optimum" step, 
(indicated by asterisks in Table II). In these cases, the initial step used was 
exactly half the stepsize indicated. 
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Table II 


Fixed Order - Variable Step 


T, = .5 x 10' 8 T 2 = .5 x 10" 13 b = .1 x 10' 10 


a 

e 

P 

Type of 
Step Mod. 

a 

Stepsize 

(Range) 

No. of 
Der. Eval. 

Est. 

Error 

6.7 

.003 

7 

H/D 


1.7 

2377 

.2 x 10' 




OPT 

.1 X 10" 9 

5.4 

744 

.6 x 10" 




OPT 

.1 x 10' 10 

3.0 

1320 

.1 x 10' 



9 

H/D 


6.7 

590 

.2 x 10' 




H/D* 


14.0 

281 

.1 x 10' 1 




OPT 

.1 x 10' 9 

12.4 

322 

.5 x 10" 




OPT - 

.1 x 10' 10 

8.4 

471 

.1 x 10' 



11 

H/D 

| 

13.4 

291 

.2 x 10" 

! 



H/D* 


23.0 

167 

.6 x 10' 




OPT 

.1 X 10' 9 

20.1 

196 

.1 x 10' 




OPT 

.1 x 10' 10 

15.1 

260 

.7 x 10" 



13 

H/D 


13.4 - 26.9 

174 

.6 x 10' 




H/D* 


23.0 

165 

.5 x 10' 




OPT 

.1 x 10' 9 

19.8 - 30.0 

155 

.8 x 10“ 




OPT 

.1 x 10' 10 

13.9 - 26.8 

159 

.7 x 10' 

1.15 

.075 

7 

H/D 


.42 

9516 

.1 x 10" 




OPT 

.1 x 10 9 

.42 

9516 

.1 x 10" 




OPT 

.1 x 10' 10 

.42 

9516 

.1 x 10" 



9 

H/D 


.42 

9514 

.2 x 10" 




OPT 

.1 x 10" 9 

.61 

6514 

.6 xl0‘ 




OPT 

.1 x 10' 10 

.42 

9490 

.2x10' 



i 

11 

H/D 


.84 

4781 

.1 x 10' 




OPT 

.1 x 10' 9 

1.04 

3849 

.8 x 10" 




OPT 

X 

o 

1 

o 

.78 

5131 

.4 x 10' 


,/D “ "Halving-doubling” 

,PT “ "Optimum Step Computation” 



Table II (Continued) 


a 

e 

P 

Type of 
Step Mod. 

a 

Stepsize 

(Range) 

No. of 
Der. Eval. 

Est. 

Error 



13 

H/D 


.84-1.68 

4257 

.4 x IO" 6 




OPT 

.1 x 10' 9 

1.20 

3314 

.8 x 10~ 9 




OPT 

.1 x nr 10 

.96 

4171 

.5 x 10' 9 

8.5 

.87 

7 

H/D 


.2 - 3.4 

3180 

.7 x 10' 8 




OPT 

.1 x 10" 9 

.2 - 9.5 

2415 

.2 X 10“ 7 




OPT 

.i x i<r 10 

.1 - 6.9 

3131 

.3 x 10' 8 



9 

H/D 


.2 - 6.7 

1374 

.7 x 10' 7 




OPT 

.1 x nr 9 

.2 - 19.0 

1137 

.7 x 10~ 7 




OPT 

.1 x 10' 10 

.2 - 14.6 

1331 

.1 x IO -8 



11 

H/D 


.2 - 13.4 

875 

.6 x io~ 7 




OPT 

.1 x 10" 9 

.2 - 29.1 

788 

.2 x 10" 7 




OPT 

.1 x io' 10 

.2 - 14.0 

907 

.5 x 10“ 8 



13 

H/D 


.2 - 26.9 

710 

.1 x 10~ 6 




OPT 

.1 x 10~ 9 

.2 - 22.4 

661 

.3 x 10" 7 


1 


OPT 

.1 X IO" 10 

.1 - 25.6 

775 

.1 x IO" 7 


(iii) For the "high e" case, (e = .87), significant gains in efficiency due to 
the larger stepsizes was effected. (No loss in efficiency due to stepsize modifi- 
cation was noted. See the results concerning the degree of control below - 
Table II.) The apparent loss in accuracy, (especially for the higher orders), 
was expected since the error control bounded the error from below as well as 
above, whereas for the corresponding results in Table I, the local error was 
insignificant for most of the integration. 

(iv) In all cases, the optimum step computation type interval modification 
yielded fewer derivative evaluations than halving-doubling, with little or no 
loss of accuracy. The effect of the magnitude of the parameter on the integra- 
tion was as expected, viz., the small value yielded an increase in accuracy at 
a cost of smaller stepsizes. 

(v) In all cases, the advantage of automatic selection of a stepsize re- 
quired to satisfy the local error criterion for the given order is evident. Also 
note that the "best" order, (yielding the fewest derivative evaluations) to use 
with the variable step control may not be the same as the "best" order to use 
for a fixed order, fixed step integration. 
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We now consider the effects of a variable-order, variable step control on 
our test orbits. The results of these integrations are given in Table III. By 
the definition of this control, the order was allowed to vary between the limits 
Lj andL 2 before any step modification was effected. In all cases, the "smallest" 
order satisfying (12) in any given stepsize was used. All integrations we per- 
formed with an initial stepsize of .42 mins, and an initial order of Lj + (L 2 ~Lj)/ 2, 
(forcing the initial order p to be in the interval (L p L 2 )). We have indicated the 
stepsize and order range which occurred during the integration as a result of 
the error control. 

Table III 

Variable Order - Variable Step 


Tj = .5 x 10' 8 T 2 = .5 x 10 1 3 8=. 1x10’ 10 a = .lxlO' 9 


a 

e 

l ,- l 2 

Type of 
Step Mod. 

Stepsize 

(Range) 

Order 

(Range) 

No. of 
Der. Eval. 

Est. 

Error 

6.7 

.003 

7-13 

H/D 

3.3 

8 

1184 

.9 x l(T 13 




OPT 

8.8 

10 

454 

.5 x 10“ 1 2 



9-13 

H/D 

6.7 

9 

588 

.2 x 10” 1 1 




OPT 

18.4 

11 

217 

.5 x 10' 10 



11 - 15 

H/D 

13.4 

11 

289 

.2 x 10' 1 1 




OPT 

24.6 

11 - 15 

175 

.2 x 10~ 8 

1.15 

.075 

7-13 

H/D 

.42 

8-10 

9513 

.9 x 10“ 9 




OPT 

.84 

10 

4767 

.5 x 10' 8 



9-13 

H/D 

.42 - .84 

9-10 

4818 

.6 x 10" 8 




OPT 

1.04 

11 

3849 

.8 x 10" 8 



11 - 15 

H/D 

.84 

1 

12 - 13 

4751 

.1 x 10' 9 




OPT 

1.2 

13 

3314 

.8 x 10' 9 

8.5 

.87 

7-13 

H/D 

.2 - 3.3 

7-13 

2744 

X 

(— 1 
o 

1 

oo 




OPT 

.2-10.8 

7-13 

2269 

.1 x 10' 7 



9-13 

H/D 

.2-13.4 

9-13 

1002 

.6 x 10' 7 




OPT 

.2-13.2 

9-13 

1067 

.3 x 10' 7 



11 - 15 

H/D 

.1-26.8 

11 - 15 

665 

.3 x 10' 6 




OPT 

.1-20.0 

11 - 15 

742 

.2 x io' 6 
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Examining these results and comparing with the previous tables we note 
the following: 

(i) For the "low e" cases, the stepsizes attained for the various orders 
selected were comparable to the stepsizes obtained for the corresponding 
orders in Table II. 

(ii) In all cases, (especially e = .87), allowing the order to vary before 
changing the stepsize resulted in a smaller "mean" stepsize and thus a larger 
number of derivative evaluations; so that if the "best" order to use with a 
variable step control is known, this control results in a more efficient integra- 
tion than the variable order-variable step control. 

(iii) In all cases, the advantage of automatic selection of both stepsize and 
order required to satisfy the given local error criterion is evident. 

In Table IV below, the effects of a variable order control alone or the aver- 
age number of predictor-corrector iterations, (and hence the total number of 
derivative evaluations), and the total error is demonstrated. The order or order 
range obtained as a result of the control is indicated. 


Table IV 

Variable Order-Fixed Step and Fixed Order-Fixed Step 


T x = .5 x 1CT 8 T 2 = .5 x 10" 13 S = .1 x 10” 10 


a 

e 

Mode 

Stepsize 

Order 

No. of 
Steps 

Avg. no. of 
p-c Iterations 

Est. 

Error 

6.7 

.003 

Vary Order 

25.0 

11 

154 

1.00 

.1 x 10” 8 



Fixed Order 

25.0 

7 

158 

1.80 

.4 x 10” 5 



Fixed Order 

25.0 

9 

156 

1.00 

.2 x 10" 6 

1.15 

.075 

Vary Order 

1.2 

11 

3327 

1.00 

.5 x 10” 7 



Fixed Order 

1.2 

7 

3331 

1.94 

.1 x 10" 4 



Fixed Order 

1.2 

9 

3329 

1.28 

.2 x 10” 5 

8.5 

.87 

Vary Order 

0.5 

5-10 

8000 

1.00 

.2 x 10” 7 



Fixed Order 

0.5 

7 

7998 

1.01 

.5 x 10” 5 
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Examining these results, we note the advantages of automatic selection of 
order both from the predictor-corrector point of view and from the accumulated 
error which results by controlling the local error. 

Finally, we wish to consider the effects of varying the "allowable range" 

(Tj, T 2 ) in the local error control. In general, the smaller the interval (Tj, T 2 ) 
is made, the greater the frequency of interval and/or order modification, and 
one would expect, particularly in the case of a variable step control, that ef- 
ficiency gains due to the larger stepsizes resulting from the control could be 
offset by the cost in changing the stepsize, should such changes occur too fre- 
quently during the integration. 

In Table V, the results of integrating one of our test orbits with a fixed 
Tj and various values for T 2 (approaching Tj as a limiting value), are tabulated. 
A variable step control (halving-doubling) was used, and the variations in the 
stepsize, along with the number of step changes and computation time (in min.) 


Table V 

Variable Step Control - Range Modification 


T, = .5 x 10' 8 S = .1 x 10‘ 10 


a 

e 

t 2 

Stepsize 

(Range) 

No. of 
Der. Eval. 

Est. 

Eri-or 

Computation 
Time (Mins.) 

8.5 

.87 

.5 x 10' 15 

.42 - 6.7 

1810 

.7 x 10" 7 

.094 



.5 x 10' 13 

.42 - 6.7 

1309 

.7 x 10" 7 

.073 



.5 x 10' 1 1 

.42 - 13.4 

886 

.6 x 10" 7 

.059 



.1 x 10" 10 

.42 - 13.4 

849 

.9 x 10" 7 

.058 



.2 x 10" 10 

.42 - 13.4 

832 

.1 x 10" 6 

.260 



.25 x 10" 10 

.42 - 26.9 

856 

.1 x 10" 6 

.810 



.3 x 10" 10 * 

.42 - 26.9 

892 

.1 x 10" 6 

2.10* 


are presented for each integration. We have indicated by an asterisk that 
integration in which a stepsize modification occurred approximately at each 
step, rendering the error control completely ineffective, since the entire 
computation time was governed by the backpoint computation. 

Examining these results, we see that as T 2 approaches Tj, gains in efficiency 
due to the larger stepsizes is first obtained, but as the range (Tj, T 2 ) becomes 
smaller, these gains are offset by the cost in changing the step, and finally, as 
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expected, causes the control to be completely ineffective. We remark, however, 
that the equations of motion used to obtain our results (see Appendix B) were not 
as complicated or lengthy as one might encounter in actual practice and in a 
more realistic situation, a much smaller interval (T 1} T 2 ) may be possible be- 
fore the computing time is completely governed by the frequency of step 
modification. 

As a final note, we remark that although we have considered in our numeri- 
cal results only a particular model of the perturbation function P it is expected 
that variations in this function, such as the inclusion of higher order gravitational 
effects, should not affect the general qualitative behavior of the local error and 
hence our results. 


SUMMARY AND CONCLUSIONS 

The problem of effecting an efficient and accurate orbital integration by a 
multistep process using the concept of controlling the local error during the 
integration has been studied. It has been shown that during the integration, the 
parameters p and h can be used to control the local error in such a way that the 
efficiency of the process is improved with no effective sacrifice in accuracy. 

In particular, it has been shown that if a sufficiently large range (T,, T 2 ) in 
the local error is allowed, and an order p is given, a variable step control can 
yield a good approximation of the optimal initial stepsize (with respect to the 
given order) automatically, and can significantly improve the efficiency of the 
process by varying this step during the integration. Moreover, if a good ap- 
proximation of the "best" order to use with a variable step - fixed order con- 
trol is known, this control effects the most efficient integration with respect to 
the various controls considered. On the other hand, a variable-order variable 
step control has the advantage that both the order and step required to satisfy 
the given local error criterion are automatically determined, and also effects 
a reasonably efficient integration. Finally, it has been seen that even a variable- 
order fixed-step control, although ineffective in minimizing the number of inte- 
gration steps, can effect a more efficient integration than no control at all, by 
minimizing the number of predictor-corrector iterations. 

We conclude by remarking that we have considered only a "local" optimiza- 
tion problem in the sense that the optimal stepsize is defined on a step-by-step 
basis as being the largest stepsize satisfying a given local error criterion at 
any given point. A more significant problem would be the consideration of the 
optimal stepsize (and/or order) distribution over the entire range of integra- 
tion, defined on the basis of a criterion on the accumulated truncation error. 
Under some restrictive conditions, (among others, that we have only a single 
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differential equation), the problem of optimizing the mesh distribution is con- 
sidered in a paper by D. Morrison.* One of the basic problems one would en- 
counter in applying such techniques in our context is the requirement that a 
"memory" of equally spaced points be available at each step during the integra- 
tion. One solution that seems worth considering is the integration of divided- 
difference interpolation polynomials, (which do not require equally spaced 
points), for use in numerical integration. 
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APPENDIX A 


DERIVATION OF THE "SUMMED" FORM OF INTEGRATION FORMULAE 


As remarked earlier, by integrating Newton's backward difference inter- 
polation polynomial, the following multi-step formulas for the numerical inte- 
gration of (1) can be derived*, (we include formulas for the velocity in the event 
that P contains x): 


(A) Equations for x - x 


Predictor: 

(Al) 

x “ 2x + x _ , 

m + 1 m ml 

= h 2 

Cr X 
0 m 



cr 0 V 2 x + cr.V 3 x + 

2 m 3 m 

• • • • + cr V k x 

k i 

Corrector: 

(A2) 

x - 2x + x , 

m + 1 m m - 1 

= h 2 

[ a *o x 


+ ^ V2 *„♦! + ••• V “ *.•>.] 


(B) Equations for x - x 

Predictor: (Bl) * m + 1 - * n = h [y 0 x m + 7l Vx m + •• + 7 k ^ k x m ] 

Corrector: (B2) x m + 1 - x m = h [y* x m + 1 + y\ Vx m + 1 + • • + 7* k vk * m + i] 

where the coefficients a , a*, y , y* are given by the following recurrence 
relationships: 


*Henrici, 1962, pp 192-195 and 291-293- 
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letting 


1 1 1 

h ra " l + 2’ + 3 + "'"' + m 




a 0 - 1, 

^0 

11 

O 

II 

l. 

* 

= i 



2 


2 



2 . 

(1) 

a 

m 

~ 1 "T h 2 

Cr , 

m l 

- T h 3 a m -2 " 


” m 

+ 2 h m+l 

(2) 

* 

- 2-k 

* „ 

2. * 


2 

i * 

a 

m 

" " 3 h 2 

m “ 1 

■ 4 h 3 °™-2 

i 

m + 2 

h 4 . i & n 
m + 1 0 

(3) 


1 

l 


l 



7 

' m 

■i'TVj 

" 3 

^m-2 ' m 

+ l 


(4) 

* 

i * 

1 

* 

l 

* 


7 

' m 

2 ^m-l 

3 

y m - 2 m 

+ l 

' 



m = 1, 2, 3, k 


v2 * m = v ( v * m ) = v (* m - K-i) = * m ' 2 " m -i + x .-i 

v k k m = v(v k -i k m ) 


It has been established (Henrici, 1960, 1962) that an algebraic equivalent 
of equations (A) & (B) , known as the "summed" form of the integration formulae, 
considerably reduces the propagation of round-off error. Formally, one can 
obtain this modification of equations (A) and (B) by applying the inverse differ- 
ence operators V -1 , V" 2 defined by V^V =1, V" 2 V 2 =1,(1 the identity), to 
both sides of these equations. In particular, by applying V 2 to both sides of 
equations (A) we obtain 


Predictor: x = h 2 fcr n V" 2 x + o-,V -1 x + cr,x + 

m+1 0 m 1 m * m 

CT-Vx + ••• + Cr. V k-2 X 

3 m km 
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Corrector: x 


m + 1 


h 2 [or; V- 2 i m + 1 + o-;V- 1 k m + 1 + + 


CT 3 V *n, + l + 


+ CT k vk " 2 * m+ l] 


If we define r S , 11 S by 

m m ^ 


V 1 * = I s m , 

m m 

v('sj = *. 


= V_1 ( ls .) = "s. 

v n s„ = ^s.) = ‘s.-'s.. 


1 W X 
1 m 


we have 


Predictor: (Al)' x m + 1 = h 2 \a Q n s + CT i S + a ,x + •• + a. V k ' 2 i 

m l m A m k 


and since 


V 2 “S.., = = ' s ... - ' s . 


X 


m + 1 


we have *S .. = x .. + *S , 

m + 1 m + 1 m 


and since 


VS-. = "S.ti - ns m = Vl 


m + 1 


we have 


“S-.. = "S. + \ tl 


m + 1 


= n S + *S + x , . 

m m m + 1 


so we have 
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Corrector: x +1 = h 2 [a* ( n S + *S + x .. ) + a* (x +1 + X S ) 

m + l 0 \ m m m + 1 / 1 \ m + 1 m / 


+ cr*x + cr^Vx + • • + Cr *yk-2 X 

°2 X m+l ^3 m + l V m + l 


or 


(A2) ' 


m + l 


h 2 [a* n S + (a* + a!) J S 

_ 0 m \ 0 1 / nr 


(cr* + a* + Cr* ) X + Cr*Vx + 
\0 1 2 j m + l 3 m + l 


+ o-*V k ~ 2 x 

k m + 


Likewise for equations (B) , applying the operator V -1 to both sides, we obtain 


Predictor: x h y_ V -1 x + y, x + y„Vx + •• + y V k_1 x 

m + l 7 0 m+ / 1 m / 2 m / k m 


Corrector: x +1 - h |~y*V X x + y*x , , + y*Vx ,-+••+ y?V k *x t 

m + l 0 m + l 7 1 m + l / 2 m+l / k m + l 


and using the definitions of J S as above, we get 


Predictor: (Bl)' x = h 

m +. 1 


y n *S + y,x + y Vx + ••• + y.V k ” 1 x 

/ 0m / lm 7 2 m 7 k n 


and as before, since V *S . = r S ^ - X S = x 

m+ 1 m + l m m + l 

we have *S .. = x .. + J S 

m+l m + l m 

and hence 


Corrector: (B2)’ x 

' m + l 


h [7l (*« + l + \) + + •• + A Vk ' lx n 


m + l 


h [rl \ + K + y \) *.+i + ' 


+ y*V k-1 x 


These are the "summed" forms of the integration formulae. (For details con- 
cerning the computational usage of these formulae, see Reference 9.) 
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We remark again that equations (Al) , (A2) and (Al)', (A2)' are algebraically 
equivalent, so that for any fixed k, any solution of (Al) , (A2) is a solution of 
(Al)', (A2)', and vica versa; and in particular, the local truncation errors 
associated with these formulas are the same; (likewise for (Bl) , (B2) and (Bl)', 
(B2) ') . 
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APPENDIX B 


EQUATIONS OF MOTION 


The equations of motion used to obtain the numerical results are given by 


where x = (x, y, z), R - ||x|| - (x 2 + y 2 + z 2 ) 1/2 , y is a constant, Fj the per- 
turbation due to the non-sphericity of the earth, and F 2 the perturbation due to 
drag. Fj = (f 1x , F ly , Fj z ) is given by 


F 


lx 


F 


ly 


F 


1 z 


jjR 2 (5s - 1) + Hz(7s - 3) + -g- (-63s 2 + 42s - 3)] 

-p jjR 2 (5s- 1) + Hz(7s ~ 3)+^ (-63s 2 + 42s - 3)J 
^ [jR 2 (5s- l) + % ( - 63s 2 + 70s - 15)J 

+ — ^ (35s 2 - 30s + 3 ) 
5R 


where s = (z/R) 2 , J, H, K are the 2 nd , 3 rd and 4 th harmonics of the earth's 
potential field. F 2 is of the form 


c d a 

2M P 


V 


V 
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where V r is the relative velocity of the satellite, p the atmospheric density at 
the satellite position; and A, M and C D are the cross sectional area, mass and 
drag coefficient of the satellite respectively. The actual numerical values of 
these constants which were used to obtain the results can be found in Refer- 
ence 9. 
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